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We present a numerically efficient technique to evaluate the Green’s function for extended two 
dimensional systems without relying on periodic boundary conditions. Different regions of interest, 
or ‘patches’, are connected using self energy terms which encode the information of the extended 
parts of the system. The calculation scheme uses a combination of analytic expressions for the 
Green’s function of infinite pristine systems and an adaptive recursive Green’s function technique 
for the patches. The method allows for an efficient calculation of both local electronic and transport 
properties, as well as the inclusion of multiple probes in arbitrary geometries embedded in extended 
samples. We apply the Patched Green’s function method to evaluate the local densities of states 
and transmission properties of graphene systems with two kinds of deviations from the pristine 
structure: bubbles and perforations with characteristic dimensions of the order of 10-25 nm, i.e. 
including hundreds of thousands of atoms. The strain field induced by a bubble is treated beyond an 
effective Dirac model, and we demonstrate the existence of both Friedel-type oscillations arising from 
the edges of the bubble, as well as pseudo-Landau levels related to the pseudomagnetic field induced 
by the nonuniform strain. Secondly, we compute the transport properties of a large perforation with 
atomic positions extracted from a TEM image, and show that current vortices may form near the 
zigzag segments of the perforation. 


I. INTRODUCTION 

Following the isolation of graphene a general class of 
two dimensional materials with widely diverse and unique 
electrical, mechanical and optical properties has been 
realized.Two dimensional materials are almost en¬ 
tirely surface and are therefore very susceptible to exter¬ 
nal influences like direct patterning^, adsorbate atoms^, 
strain^, etc. This variety of ways to alter and control 
the material properties opens a huge range of engineer¬ 
ing possibilities.^ In this context, it becomes important 
to investigate large scale disorder or patterning in re¬ 
lation to the electronic properties of graphene and re¬ 
lated two dimensional materials. From a theoretical per¬ 
spective several methods are available.^ Typically, the 
electronic structure of the system is described with a 
tight-binding type Hamiltonian and a popular approach 
is then to construct the entire system in a piece-wise 
manner using recursive Green’s functions (RGFs).^ In 
this way, we can extract the necessary terms for calcu¬ 
lating physical quantities of interest. The RGF method 
is best-suited for systems which are either finite or peri¬ 
odic in one dimension. It is frequently used for modeling 
transport, where self energies calculated using recursive 
techniques are used to attach semi-infinite pristine ieads 
to either side of a finite device region.^ Aiternativeiy, 
an efficient approach to large disordered systems is the 
real space Kubo-Greenwood approach.However, this 
method cannot include open boundary conditions and 
can only obtain average system quantities, as opposed to 
local electronic and transport properties. 

In the most common formulation, the RGF method 
treats (quasi) one dimensional systems with only two 
leads. Although variants of the method can be used for 


arbitrary geometries and multiple leads,the method 
remains limited to finite-width or periodic systems. Con¬ 
sequently, it cannot describe local and non-periodic per¬ 
turbations, or point-like probes similar to those consid¬ 
ered experimentally.An extension of recursive tech¬ 
niques, to allow efficient treatment of local properties 
in systems without periodicity or finite sizes, would al¬ 
low for easier theoretical investigation of systems which 
are computationally very expensive, or completely out of 
reach, using existing methods. 

In this paper, we develop a Green’s function (GF) 
method which is able to efficiently treat large and fi¬ 
nite sized ’patches’ embedded in an extended system, as 
shown in Fig. 1. The method combines an analytical for¬ 
mulation of the Green’s functions describing a pristine 
system^^’^^ with an adaptive recursive Green’s function 
method to described the patches. It allows for calculation 
of both local electronic and transport properties and for 
the inclusion of multiple leads and arbitrary geometries 
embedded within an extended sample. 

This patched Green’s function method exploits an ef¬ 
ficient calculation of the GF for an infinite pristine sys¬ 
tem using complex contour techniques. Using this GF, 
an open boundary self energy term can be included in 
the device Hamiltonian to describe its connection to an 
extended sample. The device region itself, containing 
nanostructures and/or leads, is then treated with an 
adaptive recursive method. We demonstrate the formu¬ 
lation using graphene, but it is generally applicable to 
all (quasi) 2D structures where the Green’s function for 
the infinite pristine system can be determined. Conse¬ 
quently, the patched Green’s function method is a versa¬ 
tile tool for efficient investigation of non-periodic nanos¬ 
tructures in extended two dimensional systems. 
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FIG. 1. The left panel shows a schematic of a computational setup containing a finite device ‘patch’ , described by Hd, 
embedded within an extended system described by the self energy X )b • The right panel shows a computational setup containing 
several device ‘patches’ of interest. 


The rest of the paper is outlined as follows: the gen¬ 
eral formalism is developed in Section IIA by calculat¬ 
ing the open boundary self energy from the pristine GF. 
In Section IIB we use graphene as an example to show 
the calculation of the pristine GF, while Section IIC dis¬ 
cusses the adaptive recursive method used to treat the 
device when including the boundary self energy. In Sec¬ 
tion III we use the developed method to study the local 
density of states of a graphene sample under the influ¬ 
ence of a local strain field. As a result, we can compare 
local density of state (LDOS) maps with pseudomagnetic 
field distributions. In this way, we show the existence of 
Friedel-type oscillations along with pseudomagnetic field 
effects in the LDOS. Finally, in Section IV, we use the 
patched Green’s function technique to demonstrate the 
existence of vortex like current patterns in the presence 
of a perforation within an extended graphene sheet. 


II. METHOD 

We consider the computational setup schematically 
shown in the left panel of Fig. 1, where a device region is 
embedded within an extended two dimensional system. 
This setup ensures that we are not including edge effects 
due to the finite-size of the simulation domain. The 
device region is described by a Hamiltonian, iT, which 
may include disorder, deformations, mean field terms or 
leads etc. This device region is embedded into an ex¬ 
tended system by applying a self energy term, To 

consider the setup in Fig. 1, we need two things: first, 
we need to construct to describe the extended part 
of the system and secondly, we need an efficient way to 
describe the device region while taking into account. 
Furthermore, the treatment of the device should be able 
to consider arbitrary geometries, including mutually dis¬ 
connected patches within the extended system, as shown 


in the right panel of Fig. 1. 

We describe the method in three steps: 

A: Derivation of the boundary self energy term, 
in terms of the pristine lattice GFs. 

B: Calculation of the real-space GF needed in the self 
energy calculation. We use graphene as an exam¬ 
ple. 

C: Implementation of an adaptive RGF method to 
build the device region(s) efficiently while includ¬ 
ing the self energy term(s) 


A. Boundary self energy 

To construct the boundary self energy describing 
the extended region in Fig. 1, we consider the simple 
graphene example in Fig. 2a. Here a central device re¬ 
gion, indicated by the dashed square, is embedded into 
an extended sheet. In this example both the extended 
area and the device region are assumed to be graphene- 
based, but the following arguments are general to any 
two dimensional material. We consider a division of the 
system into two parts: sites in the device (D) or sites 
in the extended sheet region. Furthermore, we subdivide 
the extended sheet into boundary sites (B) which are 
indicated by blue in Fig. 2 and have a non-zero Hamil¬ 
tonian element coupling them to the device region, or 
‘sheet’ sites which do not couple to the device region. 
Within a nearest-neighbour tight-binding Hamiltonian, 
the boundary sites in Fig. 2a are shown by blue symbols 
and have non-zero couplings to the device sites indicated 
by red symbols. We can now write the Hamiltonian for 
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FIG. 2. a) Shows the desired device region, indicated by the dashed square, embedded within an extended system. Red symbols 
are the edge of the device and blue symbols indicate sites in the surrounding sheet that couples to the device. We obtain the 
disconnected system discussed in the text by removing the couplings that cross the dashed line, b) Shows the corresponding 
pristine system. Again the disconnected system is obtained by removing couplings along the dashed line, c) Illustrates how 
the effect of the extended sheet on the device region is taken into account by the self energy, see Eq. (4). 


the entire system, in block matrix form, as 


Vd,b 0 \ 

Hb,B sheet I ’ (1) 

^heet,B -f^sheet / 

where the light shaded part of Eq. (1) represent an infi¬ 
nite Hamiltonian. The connections between device and 
sheet, (i.e. between the red and blue symbol sites) are 
contained in the off-diagonal blocks Vd^b and Vb,d- 
We aim to replace the infinite Hamiltonian H with a 
finite effective Hamiltonian, iTeff = Hd,d + which 
takes into account the extended sheet using a self en¬ 
ergy term To do this, we consider the connected 

system in panel a) of Fig. 2, and a disconnected system 
formed by removing the Hamiltonian elements Vd,b and 
Vb,d^ corresponding to removing couplings crossing the 
dashed line in Fig. 2a. The GFs of the connected 
and disconnected systems can be related via the 

Dyson equation, and in particular we can write the GF 
of the connected device region as 


[Hd,d 
^ — I Vb,B> 

V 0 


G (con) _ ^(dis) ^(dis)i^ ^(con) 
B,B> ~ ^D,D w ^D,B^b,D * 


( 2 ) 


Applying the Dyson equation again to obtain and 

inserting this into Eq. (2) allows us to simplify. 


= {E1-Hd,d-^b) \ (3) 

where the self energy term is given by 


‘frame’ that remains when the device is removed from the 
full system. We take advantage of this to temporarily re¬ 
place the device with a corresponding pristine region of 
the same size, as shown in panel b) of Fig. 2. The self¬ 
energy required to incorporate the finite pristine region 
into an infinite, pristine sheet is the same self energy, 5]^, 
that is required in Eq. (3). We can therefore write the 
required GF matrix, in terms of the GF of the 

infinite pristine sheet, G*^^\ These are related using the 
Dyson equation with a perturbation — Vd,B 5 

= (l + ^b]d^d,b) G^b]b- (5) 

The advantage of writing the self-energy in terms of the 
pristine sheet GFs, G^^\ and becomes clear in the 

next section, where we demonstrate an efficient method 
to calculate these two terms. It is worth noting that 
only needs to be calculated for the sites in D which 
connect to sites in B. These sites are indicated by red in 
Fig. 2 and are where the self-energy terms need to be 
added, as shown in panel c). In this way, the computa¬ 
tions only involve matrices corresponding to the edge of 
the device and not the size of the full device region as 
straight forward inversion would require. 

The calculation scheme can be summarized as follows: 

1: Calculate G^^^ and G^^j^ using the methods out¬ 
lined in Section HB. 


Sb = VD,BGf;'^VB,D. (4) 

We note that the self energy in Eq. (4) is independent 
of the considered device and depends only on GF ma¬ 
trix elements connecting sites in the pristine surrounding 


2: Calculate from Eq. (4) and Eq. (5). 

3: The finite GF for the device region, G^°^\ is given 
by Eq. (3) and can be treated using an adaptive 
RGF method, see Section HC. 
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We note that this approach does not require a spe¬ 
cific geometric shape of the device, nor does the device 
region need to be contiguous. We can treat different 
non-connected patches in an extended system, as shown 
in the right panel of Fig. 1, by extending the set D to 
include sites inside each patch and similarly expanding 
B to include sites at the boundary of each patch. The 
method presented in this section is applicable to any sys¬ 
tem where the connected, pristine GFs are easily obtain¬ 
able as demonstrated in the next section using a tight- 
binding description of graphene as an example. 

B. Real space graphene Green’s function 


on the B sublattice and Nij{z, k) = when i is on 

B and j on A. 

To simplify the notation we introduce the dimension¬ 
less k-vectors kA = 3kyao/2 and kz = V3kxaol2 such 
that f{kA, kz) = 1 + 2 cos (^kz)e^^^ , and write the sepa¬ 
ration vector in terms of the lattice vectors r = rj —r^ = 
mai + na 2 . Inserting this into Eq. (8) gives 

G°{z,r) = + J dkA J dkzNij{z,kA,kz) 

QikAi'rn-\-n)-\-ikzirn—n) 

Z _ 

z‘^ — t^(l + 4cos^(/cz) + 4 cos(/ca) cos{kz)) 

(9) 


We now turn to the calculation of the real space GF 
of the pristine system, which is needed to calculate the 
self energy, in Eq. (4) and Eq. (5). The approach 
required to calculate this quantity is demonstrated be¬ 
low for the case of a graphene sheet described with a 
nearest-neighbor tight-binding Hamiltonian, but is eas¬ 
ily generalized for other cases. 

This Hamiltonian is given by 

H=J2 ( 6 ) 

<hi> 


Eq. (9) can be solved using a two-dimensional numerical 
integration, but as we require Eq. (9) for each Green’s 
function element individually, we wish to increase the 
performance by doing one integration analytically using 
complex contour techniques. 

Eollowing the approach of Ref. 16, we use kA as com¬ 
plex variable and consider the poles, g, of the denomina¬ 
tor 


q = cos 


"4 - 1 - 4cos2 [kz)~ 

4 cos (kz) 


( 10 ) 


where the sum < i^j > runs over all nearest neighbour 
pairs and the carbon-carbon hopping integral is t ~ —2.7 
eV. The graphene hexagonal lattice can be split into two 
triangular sublattices, which we denote A and B, and 
neighbouring sites reside on opposite sublattices to each 
other. Using Bloch functions, the Hamiltonian can be 
rewritten in reciprocal space as^ 


Hk = t 



0 T 


( 7 ) 


where the matrix form arises from sublattice indexing 
within a 2 atom unit cell and we have used the defini¬ 
tion f{k) = 1 + with the lattice vectors 

di = ao(v^, 3)/2 and a 2 = ao(—V^, 3)/2 and ao the 
carbon-carbon distance. With this definition of the unit 
vectors we have the armchair direction along the y-axis 
(and zigzag along the x-axis). 

The eigenenergies and eigenstates of the system are 
easily obtained from this form of the Hamiltonian, and 
transforming back to real space allows us to write the 
desired Green’s function between sites i and j as^^’^^ 


G^j{z) 


1 /■ 2. 


(8) 


where z = E i0+ is the energy, is the area of 

the first Brillouin zone. The position of the unit cell 
containing site i is denoted by = m^ai + nia 2 with 
rrii and rii being integers. Einally we use the definition 
Nij{z,k) = z, when i and j are on the same sublattice 
and Nij{z, k) = tf{k) if i is on the A sublattice and j is 


The sign of the pole must be selected carefully to ensure 
that it lies within the integration contour, i.e. Im(g) > 
0, for contours in the positive half plane corresponding 
to the situation m + n > 0. Gare must also be taken 
with the additional phase terms that arise for opposite 
sublattice GEs. 

Using the residue theorem and integrating over a 
rectangular Brillouin zone, kA ^ [—tt; tt] and kz G 
[—7r/2; 7r/2], we finally reduce Eq. (9) to 


G«(z,r) 


i p Nij{z, g, 

47rt^ J_K ^ cos (^kz) sin (g) 

( 11 ) 


with g given by Eq. (10). A similar expression to 
Eq. (11) can be derived when using kz as first integration 
variable.The above derivation is based upon a nearest 
neighbour model, but can be generalised to also include, 
for example, second nearest neighbour terms^^ or uniax¬ 
ial strains. 

We can now use Eq. (11) to calculate the elements 
of the required GEs, and defined in Sec¬ 

tion HA. In this way, Eq. (11) can be used to fill up the 
elements of the desired matrices one at a time. Since we 
need GE matrices of size Nb x Nb and Nb x Nd, where 
Nb and Nb are the number of sites at the edge of the 
device region and in the region B, respectively, it could 
seem very ineffective to calculate one element at a time. 
However, the total number of GE elements to be calcu¬ 
lated is greatly reduced by the symmetries of the pristine 
graphene lattice. The lattice itself is six-fold symmetric 
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cell 2 


cell 3 


cell 4 


FIG. 3. The partitioning of a small graphene sample where all sites of interest are located in cell 1. Cell 2 contains all the 
sites coupling to cell 1 but which are not themselves part of cell 1. Likewise cell 3 is the sites coupling to cell 2 and so on. 
The red sites are assigned to the current cell and the lines indicate the sites still to be assigned. The previous cell and all sites 
already added are indicated by gray and white, respectively. The recursive sweep starting at the final cell and ending in cell 
1, indicated by filled arrows, gives the GFs connecting all sites of interest. We can also employ a second recursive sweep, as 
indicated by the white arrows, to obtain local properties everywhere within the device region. 


and each of these six identical wedges is in turn mirror 
symmetric, resulting in a 12-fold degeneracy of the GFs 
indexed by site separation vectors. Additionally, many 
of the required elements in and are identical. 

For instance, the onsite and nearest neighbour GF ele¬ 
ment appear many times, but only need to be calculated 
once. Taking the device region in Fig. 2 as example we 
have Nd = Nb = 20, yielding 400 individual elements for 
a brute force calculation. Instead, using symmetries and 
duplicates, we only need to calculate 38 and 42 elements 
when determining G^°^^ and G^°^\ respectively. The 
reduction becomes more significant for larger systems, as 
we generally only need to add the GF elements corre¬ 
sponding to the longest couplings. Consequently, only 
a small percentage of the GF elements need to be cal¬ 
culated individually and their values for frequently used 
separations and energies can be stored or reused to enable 
extremely fast calculation of the required self energies. 


C. Adaptive recursion for device region 

In this section we consider the device region where 
the boundary self energy can be added at the edge. 
The full GF of the device region is given by Gd = 
(^El — H]j — X)b) 5 where we have simplified the nota¬ 
tion from Eq. (3). From this GF both transport and local 
properties can be obtained. However, for most purposes 
we do not require every element of the Green’s function 
matrix element in the device region, and so to avoid a 
time consuming full matrix inversion, recursive methods 
are often applied^’^^’^^^^^. 

This section outlines an adaptive recursion method 
which efficiently includes the boundary self energy as 
well as an arbitrary device region shape and configuration 
(and number) of leads. Alternative approaches have been 
developed to treat arbitrary shaped regions with mul¬ 
tiple leads^^’^^’^^. These so-called knitting-algorithms 
add single sites at a time. They rely on a complicated 


categorizing of sites into different intermediate updating 
blocks making the theory and implementation cumber¬ 
some. Hence, we use an approach similar to the ones in 
Refs. 21-23, and employ an adaptive partitioning of the 
Hamiltonian matrix in order to bring it into the desired 
tridiagonal form suitable for recursive methods. 

Calculating physical properties generally requires cer¬ 
tain GFs connecting a specific set of sites in the device 
region. These sites of interest, for example, could be sites 
where we want to introduce defects, or couple to probes 
for transport calculations, or measure properties like the 
local density of states. We focus first on the general 
partitioning process, and then demonstrate how it can 
be quickly modified to account for the edge self-energy 
terms. We begin by placing all these sites of interest into 
recursive cell 1, as shown by the red sites in Fig. 3. We 
emphasize that the cells in this process are not of a fixed 
size and may consist of arbitrary sites which are not nec¬ 
essarily connected. Cell 2 is determined by selecting all 
the remaining unpartitioned sites which couple directly 
to sites in cell 1 via a non-zero Hamiltonian matrix ele¬ 
ment. In the example in Fig. 3, this consists of nearest 
neighbor sites of those in cell 1, which are not themselves 
in cell 1. This process is repeated until all sites in the 
device region have been allocated a cell, and is demon¬ 
strated schematically in the panels of Fig. 3 where red 
sites indicate the current cell, and dark gray or white 
sites indicate sites added to the previous cell, or to ear¬ 
lier cells, respectively. 

With the resultant block tridiagonal Hamiltonian, we 
can now employ the usual recursive algorithm, starting 
from cell n = N, so that the final step yields the required 
GF sites in cell n = 1. These terms can then be used to 
calculate observable quantities like transmission, LDOS, 
etc. Afterwards a reverse recursive sweep from n = 1 
to n = A" can be implemented to efficiently map local 
quantities like bond currents or LDOS everywhere within 
the device region^. For completeness the full recursive 
method is summarized in Appendix A including the re¬ 
verse sweep. We emphasize that the presented method 
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is not unique to graphene systems, but can be employed 
to arbitrary tight-binding-like models. 



FIG. 4. An example of the partitioning when the cell n — 1 is 
connected to the edge, and we need to include the boundary 
self-energy, T,b- In this case, all edge sites and self energy 
terms are included in cell n. The symbols are similar to Fig. 3. 


Including the boundary self-energy 

We now return to the specific case at hand where the 
recursive method outlined above needs to be adapted 
carefully to take account of the boundary self energy. In 
general 5]^ is a non-hermitian dense matrix connecting 
all edge sites of the device region. Therefore it is essential 
to assign all edge sites to the same cell. This principle is 
shown in Fig. 4. If cell n — 1 contains sites which connect 
to an edge site, then cell n must contain not only the 
edge sites directly connecting to cell n — 1, but also all 
other edge sites, as these are connected to each other via 
5]^. In this way, the cell, n + 1, must then contain all 
the sites connecting to cell n, i.e. also connecting to the 
edge, but not included in cell n. The full cell partitioning 
algorithm, including this step, is given in Appendix A. 



min LDOS max 


FIG. 5. (a) The PMF distribution calculated using the strain 
distribution in Eq. (13), dark being negative field and light 
being positive, (b-d) Real space LDOS maps for the A sub¬ 
lattice taken at the energies Ei = 0.06|t|, E 2 = 0.089|t| and 
E 3 = 0.23|t|, corresponding to energies of the first two pseudo 
Landau levels and an energy dominated by Friedel type os¬ 
cillations, respectively. The energies and the symbols corre¬ 
spond the ones used in Fig. 6. Sublattice B is similar and is 
obtained by rotating 60°. The scale bar is 5 nm. 


III. INHOMOGENEOUS STRAIN FIELDS IN 
GRAPHENE BUBBLES 

In this section, we employ the patched Green’s function 
method to a locally strained graphene system, demon¬ 
strating how it can prove a useful tool in investigating lo¬ 
cal properties of non-periodic nanostructures in extended 
two dimensional systems. 

Strain engineering has been proposed as a method to 
manipulate the electronic, optical and magnetic prop¬ 
erties of graphene.It is based on the close rela¬ 
tion between the structural and electronic properties of 
graphene. The application of strain can lead to effects 
like bandgap formation^^, transport gaps^^ and pseudo- 
magnetic fields (PMFs).^^“^^ 

Uniaxial or isotropic strain will not produce PMFs, 
although it has been shown to shift the Dirac cone of 
graphene and induce additional features in the Raman 
signal.On the other hand, inhomogeneous strain fields 
can introduce PMFs. In this case, the altered tight 
binding hoppings mimic the role of a gauge field in 
the low energy effective Dirac model of graphene. 


For example Guinea et al. demonstrated that nearly 
homogeneous PMFs can be generated by applying tri- 
axial strain. One of the most striking consequences 
of homogeneous PMFs is the appearance of a Landau¬ 
like quantization.Scanning tunnelling spectroscopy 
on bubble-like deformations see this quantization, where 
the observed pseudo-Landau levels corresponds to PMFs 
stronger than 300 

Deformations can be induced in graphene sam¬ 
ples by different techniques like pressurizing suspended 
graphene^^’^^ or by exploiting the thermal expansion co¬ 
efficients of different substrates.^^ As a result, introduc¬ 
ing nonuniform strain distributions at the nanoscale is 
a promising way of realizing strain engineering. The 
standard theoretical approach to treat strain effects em¬ 
ploys continuum mechanics to obtain the strain field. 
Several studies improve the accuracy by replacing the 
continuum mechanics by classical molecular dynamics 
simulations.The strain field can then be coupled to 
an effective Dirac model of graphene to study the genera¬ 
tion of PMFs in various geometries. In most studies, only 
the PMF distribution is considered as opposed to experi¬ 
mentally observable quantities like local density of states. 
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The framework presented in Section II enables us to treat 
the effect of strain on the LDOS directly from a tight- 
binding Hamiltonian. Consequently, we are now able to 
describe a single bubble in an extended system without 
applying periodic boundary conditions which may intro¬ 
duce interactions between neighboring bubbles. The dual 
recursive sweep then allows for efficient calculation of lo¬ 
cal properties everywhere in the device region surround¬ 
ing a bubble, enabling us to investigate spatial varia¬ 
tions in real space LDOS maps. In this section we only 
treat one nanostructure, but the patched Green’s func¬ 
tion technique efficiently handles several spatially sep¬ 
arated nanostructures, as the separation is added very 
efficiently through the self energy term. 

To account for strain within a tight binding approach 
we modify the hopping parameters.The nearest 
neighbour hopping in Eq. (6) between site i and j is given 
by the new distance, between the sites, 

tij ( 12 ) 

where the coefficient [3 = —d In t/9 In uq ~ 3.37.^^ We 
treat the deformation problem by applying an analyti¬ 
cal displacement profile {u{x^y)^ z{x^y)) matched against 
experimental data for pressurized suspended graphene. 
Here u{x^y) and z{x,y) are the in-plane and vertical dis¬ 
placements, respectively, which are induced by the ap¬ 
plied strain. For a rotationally symmetric aperture with 
radius i?, these are given, in spherical coordinates (r, 0), 
as 

z{r,e) = ho(^l - (13a) 

M(r,6') = uo^^l - (13b) 

for r < R. Here ho is the maximal height of the bub¬ 
ble and 1^0 = 1.136^0/i? is a constant relating the out- 
of-plane and in-plane deformations.^^ We note that this 
profile gives rise to a sharp edge at r = i?, and many of 
the features we discuss below emerge from the strongly 
clamped nature of this bubble type. 

As shown in Appendix B, rotationally symmetric 
strain profiles give rise to threefold symmetric PMFs in 
the effective Dirac model. This is shown in Fig. 5a for 
the strain profile considered in Eq. (13). As discussed 
in earlier studies,we get an asymmetric sublattice 
occupancy such that the LDOS of each sublattice has 
a threefold symmetric distribution following the PMF 
while rotated 60° compared to the opposite sublattice. In 
all calculations below, we therefore show only one sub¬ 
lattice, as the result for the opposite sublattice can be 
obtained by a 60° rotation and the total pattern is a 
superposition of both.^^’^"^ 

Comparing the PMF distribution in Fig. 5a with the 
calculated LDOS maps at different energies in Fig. 5b-d 
for a bubble of radius R = 10 nm and height ho = 3 
nm, we immediately notice that the threefold symmetry 



Energy / \t\ 

FIG. 6. The LDOS as a function of energy for the three 
positions indicated in Fig. 5 and for the average of the ‘slice’ 
of the bubble region containing the symbols. The dashed lines 
indicate the LDOS without the bubble. The curves are shifted 
with respect to each other to increase visibility. 



Energy / \t\ 

FIG. 7. The difference in LDOS as a function of energy for the 
point indicated with a triangle on Fig. 5. We show both the 
full calculation (full line) and an artificial system containing 
only the perturbation for a small region at the edge of the 
bubble (dashed line). We adjust the average hopping constant 
in the calculation of the artificial system to match the full 
calculation. Inset: The peak energies 1-4 as a function of ^/n, 
where n is the peak number. 

is also present in the LDOS maps. However, the spatial 
LDOS maps have significant additional details compared 
to the PMF distribution. 

In Fig. 6 we calculate the energy dependent LDOS at 
the positions indicated by symbols (square, circle and 
triangle) in Fig. 5. We first consider the average of the 
LDOS within the ‘slice’ containing the symbols, shown 
by the bottom (red) curve in Fig. 6. Two distinct oscilla¬ 
tion types are observed, and we argue that these can be 
divided into Friedel-type and PMF-induced features. At 
high energies in particular we notice regularly space oscil¬ 
lations with an approximate period of hv^'K/2R. These 
are consistent with Friedel-type oscillations related to the 









size of the structure and emerging from interferences be¬ 
tween electrons scattered at opposite sides of the bubble. 
An exact treatment needs to take into account the renor¬ 
malized Fermi velocity, due to the average change 
in bond length.At lower energies we observe distinct 
peaks which are not equally spaced (the first two ap¬ 
pear at El and E 2 ). We will show that these are due to 
pseudomagnetic effects and we refer to them as pseudo 
Landau levels. 

Besides the Friedel oscillation associated with the bub¬ 
ble radius, we also have similar oscillations associated 
with the distances to different edges of the bubble. These 
features are highly position dependent, and explain the 
differences between the three single position curves in 
Fig. 6. When considering the average, these position 
dependent oscillations are washed out (bottom curve in 
Fig. 6), leaving only the oscillation dependent on the 
structure size. However, at individual positions these os¬ 
cillations can have a considerable impact. Returning to 
the individual position STS curves in Fig. 6, we note that 
the peak at E 2 is only dominant for the points indicated 
by the square and triangle. It is suppressed by Friedel- 
type interferences at the circle point, which is also clear 
from the LDOS map in Fig. 5c. 

The amplitude of the Friedel-type oscillations is de¬ 
termined by the strength of scattering near the bubble 
edges. The clamped edge implied by the strength profile 
in Eq. (13) gives rise to significant strain fields along this 
edge, leading to a sharp, strong perturbation. More re¬ 
alistic profiles calculated from molecular dynamics also 
indicate strong perturbations near the edges of clamped 
bubbles.Our results indicate that edge scattering ef¬ 
fects may significantly affect LDOS behavior in clamped 
bubble systems and even mask PMF-induced features. 

To treat the oscillations due to the feature size and 
edge sharpness in more detail, we calculate the LDOS for 
an artificial system only taking into account the strain 
field along a small ring around the edge, see Fig. 7 
(dashed red line). In this way, only Friedel-type fea¬ 
tures are expected within the structure. If we com¬ 
pare to the full calculation (full black line in Fig. 7), 
we notice that the oscillations at higher energies are 
present in both calculations, whereas the sharp peaks 
are only present in the full calculation. This confirms 
the Friedel nature of the higher energy oscillations and 
suggests the lower energy peaks are due to an alterna¬ 
tive mechanism. To confirm that the sharp peaks are 
due to pseudomagnetic effects, we compare the peak po¬ 
sitions to the standard for m expected f or Landau levels 
in graphene En = where eo is the 

electron charge, Bg is the magnetic field and n is the peak 
number.^ The peaks labelled 1-4 in Fig. 7 display the y/n 
dependence characteristic of Landau levels in graphene, 
as shown in the inset of Fig. 7. The size of the PMF can 
furthermore be inferred to be Bg ^ 30 T from the inset. 

To conclude, we discussed how the features in the 
LDOS spectra of clamped graphene bubbles can be ex¬ 
plained by a combination of size-dependent scattering 
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FIG. 8. (a) The transmission as a function of energy for 

a dual probe setup with an antidot in between the probes 
as schematically shown in the inset. The distance between 
the probes are 200 nm and the antidot with purely zigzag 
edges has side length R = 48ao ~ 6.8 nm. The shaded area 
corresponds to the LDOS around the edge of the antidot. 
antidot 


and PMF-induced effects like pseudo Landau quantiza¬ 
tion. Significant strain fields near the edge of the struc¬ 
ture give rise to strong Friedel-type oscillations in the 
LDOS and these oscillations envelope the effect of a PMF. 
We must therefore be careful to distinguish between the 
two type of oscillations when investigating the electronic 
effects of PMFs induced by inhomogeneous strain fields. 


IV. VORTEX CURRENTS NEAR 
PERFORATIONS 


In this section we investigate local transport proper¬ 
ties near antidots {i.e. perforations) in a graphene sheet. 
Periodic arrays of antidots have been studied as a way 
to open a bandgap in graphene^^”^^ or to obtain waveg- 
uiding effects.Furthermore, a single perforation in 
a graphene sheet has been considered as a nanopore for 
DNA sensing. 

Several studies show that the electronic struc¬ 
ture of antidots is closely related to the exact 
edge geometry.^^’^^’^^ Experimental fabrication tech¬ 
niques like block copolymer^’^^’^^ or electron beam 
lithography,^^”^^ inevitably lead to disorder and imper¬ 
fect edges. However, it may be possible to control the 
edge geometry of the antidot by heat treatment^^’^^, or 
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selective etching. 

Motivated by the interest in how current flows in an- 
tidot systems, we apply the patched GF method to a 
single perforation in a graphene sheet. The method al¬ 
lows us to study the perforation with no influence from 
periodic repetition or finite sample size. Additionally, the 
combination of recursive methods and a boundary self¬ 
energy allows for investigation of antidot sizes realizable 
experimentally.In fact we consider both an exam¬ 
ple antidot with perfect edges and an exact structure 
found from high resolution transmission electron micro¬ 
scope (TEM) images using pattern recognition. 

To investigate current on the nanoscale, recent ex¬ 
periments have realized multiple STM-systems.^^’^^’^^’^^ 
These allow for individual manipulation of several STM- 
tips in order to make electrical contact to the sample near 
the considered nanostructure. Theoretically, we previ¬ 
ously considered multiple STM setups allowing for both 
fixed and scanning probes.The method presented 
here allows for not only transmission calculations but 
also calculation of local electronic and transport prop¬ 
erties in the presence of multiple point probes. At the 
same time large separations between the different probes 
and/or nanostructures are easily included as additional 
separation is achieved in a very computationally efficient 
manner through the self energy term connecting multi¬ 
ple patches. The combination of large spatial separation 
between features, while still enabling calculation of local 
electronic and transport properties, can prove a useful 
tool in investigating extended two dimensional systems 
where we take special interest in a particular region of 
the extended sample. 

In order to consider transmissions and current pat¬ 
terns, we add leads to the system through inclusion of 
a lead self-energy term, EF = V^g^{ri — where 

is the coupling element between the device site i and 
the lead. To model the structureless lead, we use the 
surface GF of a single atomic chain, as this has a con¬ 
stant DOS in the considered energy range. The distance 
dependence in g^{ri — Vj) is necessary to avoid an un¬ 
physical coupling between different lattice sites via the 
lead. We therefore add a l/\ri — rj|-dependence for the 
off-diagonal terms^^ where 7 ^ rj, as appropriate for a 
structureless three-dimensional free electron gas.^^ 

First, we consider a zigzag-edged antidot with side 
length R = 48a ^12 nm, where a is the length of the 
graphene unit cell and a = ^/3ao = 2.46 A. This is com¬ 
parable to experimental sizes where sub- 20 -nm feature 
sizes have been reported.The antidot is between 
two probes placed 200 nm apart, as shown schematically 
in the inset of Fig. 8 a. The main panel of Fig. 8 a shows 
the transmission as a function of energy for this dual 
point probe setup. We note the distinct transmission 
peaks. As explained in Ref. 56, these peaks are related 
to localized states along the zigzag edges. As a conse¬ 
quence, we notice the correspondence between the peaks 
in the transmission and the peaks in the LDOS around 
the edge, see shaded area in Fig. 8 a. 




FIG. 9. An actual perforation is obtained from high resolution 
TEM images through pattern recognition and we consider the 
vortex like current paths forming around the perforation as 
certain energies, (a) Shows the structure of the perforation as 
well as an indication of the probe position (in the actual cal¬ 
culations the probes are 200 nm apart). The indicated areas 
corresponds to the zooms in (c) and (d). (b) The transmis¬ 
sion for the dual probe setup. The shaded area indicate the 
average LDOS around the edge of the antidot. Furthermore, 
the energies I and II corresponding to the energies used on 
(c) and (d), respectively, (c-d) Bond current maps taken at 
the energies I and II, respectively, and shown at the positions 
indicated on (a). 


Next, we calculate the bond currents from the top 
lead. The bond current between site i and j from 
lead L are calculated, as explained in Appendix A, by 
JIj = —Hijlm[G^T^G'^].^/h, where Hij is the Hamilto¬ 
nian matrix element connecting site i and j. The bond 
currents around the zigzag antidot for the energies indi¬ 
cated in Fig. 8 a are shown in Figs. 8 b and 8 c. In this way, 
we see that the transmission dips are related to vortex 
like current paths. These vortex paths create a larger ‘ef¬ 
fective size’ for the antidot at this energy, characterized 
by a region around the antidot avoided by the current 
paths. On the other hand, at the transmission peaks the 
current passes near to the antidot edge. 

The antidot considered in Fig. 8 , although of realis¬ 
tic size, is an idealization, as experimental perforations 
will inevitably contain imperfections. To consider a more 
realistic case, we turn to a perforation observed in exper¬ 
imental TEM images. Using pattern recognition^^the 
positions of the individual carbon atoms are obtained 
from high resolution TEM images (see Fig. 9 of Ref. 
66 ). Pristine graphene is added around the experimen¬ 
tally obtained perforation to obtain the system shown in 
Fig. 9a. From the transmission (see Fig. 9b), we notice 
that peaks are still present, but broadened by the disor¬ 
der. Gonsidering the two energies I and II in Fig. 9b and 
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comparing their spatial current maps, we find that cer¬ 
tain positions around the antidot are responsible for the 
additional backscattering causing the transmission dips. 
Dip I corresponds to a vortex pattern at the left side of 
the antidot (see Fig. 9c), whereas the dip at II is caused 
by a vortex pattern at the bottom of the antidot (see 
Fig. 9d). This result suggests that electrons at differ¬ 
ent energies see a different effective perforation size and 
shape and are scattered accordingly. 


V. CONCLUSION 

We have expanded the standard recursive Green’s 
function method to calculate local and transport proper¬ 
ties enabling calculations in extended non-periodic sys¬ 
tems. We exploit an efficient calculation of the pris¬ 
tine two-dimensional GF using complex contour meth¬ 
ods. Once calculated, the pristine GFs are used to de¬ 
termine a boundary self energy term describing the ex¬ 
tended system. In this way, we can treat a finite device 
region embedded within an extended sample. 

We first demonstrated how this approach is able to effi¬ 
ciently treat the electronic properties of strained bubbles 
in an extended graphene sheet. Gonsidering a clamped 
bubble, we have shown that the finite size gives rise to 
Friedel-type oscillations in the density of states. This ef¬ 
fect mixes with any pseudomagnetic effects arising from 
the strain field. We show that the edge effects can cloud 
pseudomagnetic signatures in the LDOS by adding addi¬ 
tional structure which is not directly related to pseudo¬ 
magnetic effects. 

Secondly, we showed how finite leads can be added to 
a patched device region to efficiently calculate transport 
properties for spatially separated features, while still be¬ 
ing able to map local properties in various parts of the 
system. In particular, we investigated the current flow 
around perforations of a graphene lattice. Both idealized 
geometries and experimental geometries obtained from 
high resolution TEM images were considered. The trans¬ 
missions show distinct dips caused by localized states 
along zigzag segments of the perforations. The transmis¬ 
sion dips were associated with vortex-like current paths 
formed near the perforation edges. 

We have demonstrated the versatility of this novel ap¬ 
proach to the popular recursive GF method. The method 
allows for calculation of the same local and transport 
properties as standard methods, but adds the ability to 
treat large non-periodic structures embedded in extended 
samples. We can extend the present method beyond 
nearest neighbor and to relevant alloys like hBN or tran¬ 
sition metal dichalcogenides. We therefore predict that 
the patched Green’s function method will prove a valu¬ 
able tool in the investigation of nanostructures in two 
dimensional materials. 
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Appendix A: Recursive Algorithm 

To obtain a tridiagonal Hamiltonian we let cell n = 1 
contain all sites of interest. Then following the algorithm 
outlined below we assign all sites into cells. 

1: Let {n} denote all sites in cell n and 
{^^unassigned^^} denote all sites not yet assigned 
to a cell. 

2 : Eind all sites j for which H^j 7 ^ 0 where n G {n} 
and j G unassigned^^}. Denote these sites {n + 

!}• 

2a: If {n + 1} contains an edge site, then all remaining 
edge sites are added to {n + 1 }. 

3: Sites in + 1 } are removed from {^^unassigned^^} 

4: Repeat 1-3 until all sites are assigned to a cell. 

Step 2 a is included if we require an edge self energy term 
Eb as described in Section II. 

Assuming the block tridiagonal partitioning obtained 
from the algorithm above, we make an update sweep 
starting from cell n = A, as shown schematically in 
Eig. 10. The steps are calculated using the recursive 
relations^ 

9 n,n = {E — Hn,n) , (Ala) 

9n,n — i^E i 

(Alb) 

M 

gi.l = {E- ifi,i - V,292,2^2,1 - ^ 

m=l 

(Ale) 

where one of the iTn,n terms includes the self energy and 
E^ad l^iaas are included if we calculate transmission. Af¬ 
ter the sweep is complete, the fully connected GE of cell 
n = 1 is obtained as Gip = gi^i. As all sites of interest 
are placed in this cell, we can now calculate observables 
involving these sites. Eor example we calculate transmis¬ 
sion, 7 l,lg between lead L and L’ using these GEs. 

rL,L'{E) = (A2) 

where = i(S^ - S^t) and Gl,l' is the re- 

tarded (advanced) GE connecting the two leads L and 
L’. 

In order to obtain other blocks of the full GE matrix, 
we need to store the GE matrix, gn,n^ for each cell as 
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FIG. 10. Top row: recursive sweep going from cell n = N to n = 1. Light gray indicate blocks that are stored for the reversed 
sweep and dark gray indicate blocks of the full GF. The inset shows an illustration of how the different blocks correspond to 
neighboring cells in the device region. Bottom row: reversed recursive sweep going from n=l 1 to n = iV showing how this 
sweep can obtain both diagonal and off-diagonal blocks. 


we sweep from n = N to n = 1. The stored blocks are 
shown in light gray on Fig. 10. 

To obtain the LDOS at site i, pa = —Im(Gii)/7r, we 
need the diagonal of the GF matrix. We calculate the 
block diagonal from a reversed sweep from n = 1 to 
n = see Fig. 10. The reversed sweep uses the block 
diagonals, gn,n^ from the first sweep to calculate the full 
diagonal GF, G, 


where u{x,y) is the in-plane deformation field and z{x^y) 
is the out-of-plane deformation.^^ 

A general two dimensional strain field, eij{x^y)^ leads 
to a gauge field in the effective Dirac Hamiltonian of 
graphene^^’^^ 


2 eao 


I ^xx ^yy 
\ ‘^^xy 


(B2) 


^n,n — 9n,n T 9n,n^n,n —1^71—1,n—l^n — l,n9n,n' 

Finally, we want to obtain bond currents for the state 
leaving a lead L. This can be calculated by = 

—Hijlm[Gi^iTi iG\j]/h. Remembering that the leads 
are assigned to cell n = 1, we need the off-diagonal 
blocks, Gi^n and Gn,i, in order to obtain bond currents. 
Using the stored GFs from the first sweep we can calcu¬ 
late the needed off-diagonals, 

^l,n ~ ^l,n—l^n—l,n9n,n') (A4a) 

^n,l — 9n,n^n,n—lGfi—irf^. (A4b) 

Appendix B: Pseudomagnetic field for rotational 
symmetric strain field 

The strain tensor is generally given as 


which in turn gives a PMF 

Bs = V X A = dxAy - dyAx. (B3) 

Eqs. (B2) and (B3) imply that the x-axis is chosen along 
the zigzag direction of the graphene lattice. 

Now restricting ourselves to rotationally symmetric de¬ 
formations, u{r) = Ur and z{r) = z, while using polar 
coordinates (r, 0) yields 

Bs = _ drgir)] sin(36»), (B4) 

2eao \ r ) 

with ^(r) = drUr — Urjr + \ (drzY• We notice from 
Eq. (B4) that the PME for a rotationally symmetric dis¬ 
placement field is always 6-fold symmetric. On the other 
hand, the magnitude depends on both the in-plane and 
out-of-plane displacement. 

Considering the displacement field in Eq. (13) we now 
obtain a PME of the form. 


- djUi + diUj + {diz){djz) , ij = x,y, (Bl) 


Bs 


2eaoR‘^ 


sin(36>). 


(B5) 
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Taking into account the scaling uq oc h^jR^ we obtain 


a final scaling of the PMF with the size of the bubble, 
Bs oc hl/R^. 
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